import numpy as np
import load
load = reload( load )
# constants
TOL = 0.001
thetaTOL = 10.*TOL
def rec2cyl( x, y, z):
    thetarad = np.arctan2( y, x )
    r     = x / np.cos( thetarad )
    theta = np.rad2deg( thetarad )
    z     = z
    return r, theta, z

def cyl2rec( r, theta, z ):
    x = r * np.cos( np.deg2rad( theta ) )
    y = r * np.sin( np.deg2rad( theta ) )
    z = z
    return x, y, z

def find_fb_load( cczload ):
    frmlen = len(cczload)
    fb_load = 0.
    found_fb_load = False
    for i in range(frmlen):
        zload = cczload[i]
        if abs(zload) < abs(fb_load):
            found_fb_load = True 
            break
        if not found_fb_load:
            fb_load = zload
    return fb_load


class PLoad(object):
    
    def __init__( self, theta=0., pt=0.5, pltotal = 1. ):
        self.name   = ''
        self.index  = None
        self.theta  = theta
        self.pt     = pt       #ratio along cone/cylinder height
        if abs(pltotal) < 0.1 * TOL:
            pltotal = 0.1 * TOL
        self.pltotal = pltotal #resultant pload
        self.plradial = None   #component radial direction
        self.plx = None        #component x      direction
        self.ply = None        #component y      direction
        self.plz = None        #component z      direction
        self.x = None          #coordinate x
        self.y = None          #coordinate y
        self.z = None          #coordinate z
        self.vertice = None
        self.ploadconf = None
        self.poly_cir = None
        self.poly_mer = None
        self.xfit_cir = None
        self.xfit_mer = None
        self.ehr_cir = None
        self.ehr_mer = None
        self.displ   = None
        self.node = None
        self.cross_section = None
        self.meridian = None

    def get_xyz( self ):
        r, z = self.ploadconf.conecyl.r_z_from_ppoint( self.pt )
        return cyl2rec( r, self.theta, z )

    def rebuild( self, alpharad, betarad ):
        if abs(self.pltotal) < 0.1 * TOL:
            self.pltotal = 0.1 * TOL
        if self.theta > 180. and self.theta <= 360.:
            self.theta = self.theta - 360.
        self.plradial =  self.pltotal * np.cos( alpharad + betarad )
        self.plz = -self.pltotal * np.sin( alpharad + betarad )
        self.plx = -self.plradial * np.cos( np.deg2rad( self.theta ) )
        self.ply = -self.plradial * np.sin( np.deg2rad( self.theta ) )
        self.x, self.y, self.z = self.get_xyz()
        self.name = 'pt_%1.2f_PLvalue_%04d' \
                    % ( self.pt, int(10*self.pltotal))
        self.name = self.name.replace('0.','0')
   
    def get_node( self ):
        if self.node == None:
            for meridian in self.ploadconf.conecyl.meridians:
                if abs(self.theta - meridian.theta) < thetaTOL:
                    meridian.ploads.append( self )
                    self.meridian = meridian
                    break
            for cross_section in self.ploadconf.conecyl.cross_sections:
                if abs(self.z - cross_section.z) < thetaTOL:
                    cross_section.ploads.append( self )
                    self.cross_section = cross_section
                    break
            z_index = self.cross_section.nodes[1].z_index
            self.node = self.meridian.nodes[ z_index ] 
            self.node.pload = self

        if self.node == None:
            print 'WARNING - node not found'
        return self.node

    def calc_ehr_and_displ( self, step_name='',\
                            cir_criterion = 'zero_radial_displ',
                            cir_threshold = 0.0125,
                            only_cir = True ):
        if self.pltotal < TOL:
            self.ehr_cir = 0.
            self.ehr_mer = 0.
            return 
        if step_name == '':
            step_name = self.ploadconf.conecyl.step1Name
        self.get_node()
        cc = None
        if self.ploadconf <> None:
            if self.ploadconf.conecyl <> None:
                cc = self.ploadconf.conecyl
        if cc == None:
            print 'WARNING - concecyl not found for pload'
        #
        # circumferential equivalent hole radius
        # 
        nodes = self.cross_section.nodes
        nodeids = [ node.id() for node in nodes ]
        plnindex  = nodeids.index( self.node.id() )
        plndispl  = self.node.dr[ step_name ][ -1 ]
        if cir_criterion == 'zero_radial_displ':
            # finding point
            x = []
            y = []
            for node in nodes:
                x.append( node.theta - self.theta )
                y.append( node.dr[ step_name ][ -1 ] )
            thr = cir_threshold
            for i in range( len(y)-2, -1, -1 ):
                yi = y[i]
                if abs(yi) > thr:
                    xi  = x[i]
                    xi1 = x[i+1]
                    yi1 = y[i+1]
                    # linear interpolation
                    radius_theta = (thr-yi)*(xi1-xi)/(yi1-yi) + xi
                    break
        if cir_criterion == 'first_cross':
            # reading points to polyfit
            xj = []
            yj = []
            order = 6
            xj.append( nodes[plnindex].theta - self.theta )
            yj.append( nodes[plnindex].dr[ step_name ][ -1 ] )
            for j in range( plnindex + 1, len(nodeids) ):
                node_1 = nodes[ j-1 ]
                node   = nodes[ j   ]
                dri_1  = node_1.dr[ step_name ][ -1 ]
                dri    = node.dr[ step_name ][ -1 ]
                xj.append( node.theta - self.theta )
                yj.append( dri )
                if dri_1 > dri and (j-plnindex) > order:
                    break
            # creating polynom
            if len(xj)-1 < order:
                order = len(xj) - 1 
            poly = np.poly1d( np.polyfit( xj, yj, order ) )
            # curve to find theta for radial_displ = 0
            init  = int(10*min(xj))
            final = int(10*max(xj))
            xfit = np.array( [j*0.1 for j in range( init, final, 1 )] )
            roots = []
            for root in np.roots( poly ):
                if np.imag( root ) == 0:
                    root = np.real( root )
                    if  min(xj) < root\
                    and max(xj) > root: 
                        roots.append( root )
            if roots == []:
                for root in np.roots( poly.deriv() ):
                    if np.imag( root ) == 0:
                        root = np.real( root )
                        if  min(xj) < root\
                        and max(xj) > root:
                            roots.append( root )
            roots.sort() 
            radius_theta = roots[0]
            if roots[0] < cc.plyt:
                radius_theta = roots[1]
            #
            # storing values
            self.xj_cir = xj
            self.yj_cir = yj
            self.poly_cir = poly
            self.xfit_cir = xfit
        r_local, tmpz  = cc.r_z_from_ppoint( self.pt )
        try:
            rad = np.deg2rad(radius_theta - self.theta)
        except:
            print 'ERROR - Diameter could not be calculated'
            return False
        self.ehr_cir  = float( rad * r_local)
        #
        # meridional equivalent hole radius
        # 
        if only_cir:
            self.ehr_mer = 0.
            return
        nodes = self.meridian.nodes
        nodeids = [ node.id() for node in nodes ]
        plnindex = nodeids.index( self.node.id() )
        plndispl = self.node.dr[ step_name ][ -1 ]
        xj = []
        yj = []
        # reading points to polyfit
        order = 6
        xj.append( nodes[plnindex].z - self.z )
        yj.append( nodes[plnindex].dr[ step_name ][ -1 ] )
        for j in range( plnindex + 1, len(nodeids) ):
            node_1 = nodes[ j-1 ]
            node   = nodes[ j   ]
            dri_1  = node_1.dr[ step_name ][ -1 ]
            dri    = node.dr[ step_name ][ -1 ]
            xj.append( node.z - self.z )
            yj.append( dri )
            if dri_1 > dri and (j-plnindex) > order:
                break
        # creating polynom
        if len(xj)-1 < order:
            order = len(xj) - 1 
        poly = np.poly1d( np.polyfit( xj, yj, order ) )
        # curve to find theta for radial_displ = 0.1  
        init  = int(10*min(xj))
        final = int(10*max(xj))
        xfit = np.array( [j*0.1 for j in range( init, final, 1 )] )
        polyxxx = poly.deriv().deriv().deriv()
        #
        roots = []
        for root in np.roots(polyxxx):
            if np.imag( root ) == 0:
                root = np.real( root )
                if  min(xj) < root\
                and max(xj) > root:
                    roots.append( root )

        # storing values
        self.xj_mer = xj
        self.yj_mer = yj
        self.ehr_mer = float(min( roots ))
        self.poly_mer = poly
        self.xfit_mer = xfit
        self.displ    = plndispl

class PLoadConf(object):
    '''
       Configurations of pload
       - each PLoadConf may have many PLoad objects, applied at different
         regions of the cone / cylinder, and at different values
    
    '''
    def __init__( self ):
        self.ploads = []
        self.name = ''
        self.conecyl = None
    
    def add_pload( self, theta, pt, pltotal ):
        import load
        load = reload( load )
        pload = load.PLoad( theta, pt, pltotal )
        pload.ploadconf = self
        self.ploads.append( pload )

    def rebuild( self ):
        # ploads
        for i in range(len(self.ploads)):
            pload = self.ploads[i]
            pload.index = i
            pload.rebuild( self.conecyl.alpharad, self.conecyl.betarad )
        # name
        if self.ploads == []:
            self.name = 'NO_PLoads'
        elif len( self.ploads ) == 1:
            self.name = self.ploads[0].name
        elif len( self.ploads ) >  1:
            self.name = '%02d_PLoads_PLvalue_%04d' \
                        % (len(self.ploads), int(self.ploads[0].pltotal))
        else:
            print 'DEBUG',self.ploads
        if self.conecyl.imperfection_activate:
            num = self.conecyl.imperfection_scaling_factors[0]*10000
            self.name = 'imp_%05d' % num
                

class PLoadCurve(object):
    
    def __init__( self ):
        self.ploadconfs = []
        # real holes
        self.holes = []
        # will be equivalent holes with a given equivalent
        #     diameter corresponding to the buckle diameter
        self.buck_loads = []
    
